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ABSTRACT 

Context. We develop an optimization principle for computing stationary MHD equilibria. 

Aims. Our code for the self-consistent computation of the coronal magnetic fields and the coronal plasma uses non-force-free MHD equilibria. 
Previous versions of the code have been used to compute non-linear force-free coronal magnetic fields from photospheric measurements. The 
program uses photospheric vector magnetograms and coronal EUV images as input. We tested our reconstruction code with the help of a 
semi-analytic MHD-equilibrium. The quality of the reconstruction was judged by comparing the exact and reconstructed solution qualitatively 
by magnetic field-line plots and EUV-images and quantitatively by several different numerical criteria. 

Methods. Our code is able to reconstruct the semi-analytic test equilibrium with high accuracy. The stationary MHD optimization code 
developed here has about the same accuracy as its predecessor, a non-linear force-free optimization code. The computing time for MHD- 
equilibria is, however, longer than for force-free magnetic fields. We also extended a well-known class of nonlinear force-free equilibria to the 
non-force-free regime for purposes of testing the code. 

Results. We demonstrate that the code works in principle using tests with analytical equilibria, but it still needs to be applied to real data. 
Conclusions. 
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. . . 1. Introduction 

Understanding many physical phenomena in the solar corona 
requires detailed knowledge of the properties of the coro- 
nal magnetic field and plasma. Usually neither the coronal 
magnetic field nor the plasma density, pressure, temperature, 
and flow speed are directly observable. Direct observations of 
the coronal magnetic field are difficult because the high tem- 
perature broadens the line-profile orders of magnitude above 
the Zeeman splitting. The optically-thin coronal line emis- 
sion has a line-of-sight integrated character, which complicates 
the computation of any plasma quantities. A way out of this 
dilemma is the use of coronal models, which are fed with as 
much observations as possible. The model assumptions, e.g. an 
MHD model, constrains the non-observed quantities. Often it 
is not the full set of MHD equations used, but a subset. 

The simplest approach for coronal magnetic fie ld calcu- 
lations is to assum e a potential magnetic field (e.g. Schmidtl 
1964t ISemellll967t) . The only observational input required is 
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the line-of-sight photospheric magnetic field, such as observed 
by SOHO/MDI. These source surface potential field models 
provide a first impression regarding the global coronal mag- 
netic field, e.g. regarding the location of coronal holes and ac- 
tive regions. Details of the magnetic field structure are often 
not well-approximated by potential field models, particularly 
in active regions. Approaches using force-free magnetic fields 
(with electric currents parallel to the magnetic field) show re- 
sults which are significantly better. The low plasma beta in the 
corona (J3 ~ 10 4 ) justifies that approach in the low and middle 
corona, but not in the photosphere (j3 « 1) where non mag- 
netic forces like pressure gradients become important. A popu- 
lar si mplification of force-free fields are lin ear force-free fields 
(e.g. IChiu & Hiltonlll977l; ISeehafet|[l978l) ) where the electric 
current flow is parallel to the magnetic field with a global con- 
stant of proportionality a. 

A more sophisticated approach is to allow a to change 
in space, the so-called nonlinear force-free approach. The 
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calculatio n of no n- linear force-free fi elds (e .g. Sakur ailll981 
Wu et all 1 19901: I Wheatland et al J l2000t lYan & Sakurail 
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Wiegelmann et al.1 120061: 2. Basic equations 



Schriiver et alj 120061) is 



particularly challenging due to the intrinsic nonlinearity of 
the underlying mathematical problem. From an observational 
point of view the nonlinear reconstruction is also more de- 
manding because photospheric vector magnetograph data are 
required to determine the boundary conditions. 

A comparison of t he magne t ic loo ps measured by 
Solanki et all J2003b and lLagg et al.1 d2004 with diff erent ex- 
trapolated field models by I Wiegelmann et al. I d2005bl) revealed 
that linear force-free fields are better than potential fields, but 
non-linear force-free models are even more accurate. Usually 
the coronal magnetic field cannot be measured directly, also 
some progress has been made recently by using magnet ically 
sensitive coronal line emission, e.g., by iLin et al.1 (120041) who 
did spectropolarimetric measurements of the forbidden Fe XIII 
1075 nm infrared coronal line. The influence of the coronal 
magnetic field onto the emitting plasma is however observed 
routinely. 

Images of the coronal plasma emission are obtained by the 
radiation in different wave lengths, e.g. in EUV for SOHO/EIT 
and TRACE. The emission is obtained from various elements, 
e.g. Fe XII or Fe IX, with different instrument channels sen- 
sitive to emission generated at different plasma temperatures. 
The corresponding images have a line-of-sight character be- 
cause of the corona is optically thin. A rather good approx- 
imation of the coronal temperature is sometimes possible, 
because the radiation only occurs in a specific temperature 
range. Doppler shifts in the line profile observed with e.g. 
SOHO/SUMER provide insight into the plasma flow speed in 
the line-of-sight direction. The high electrical conductivity of 
the coronal plasma ensures that the plasma is frozen into the 
magnetic field. This basically allows us to see (the effects of) 
magnetic field lines in EUV-images and even use the visible 
plasma radiation to improve coronal magnetic field models, 
e.g. to sp ecify the optimal value a for linear force-free models 
(see e.g. Wiegelmann & Neukirchl l2002t ICarcedo et al 1 120031: 



Marsch et al.ll2004U Wiegelmann et alJl2005al) 



For the low B corona, it is also helpful that the back reac- 
tion of the p l asma onto the magnetic field can be neglected. 



Marsch et al.l (120041) used a linear force-free coronal magnetic 



field model and Doppler maps from SO HO/SUMER to in 



vestiga te the plasma flow in ac tive regions. IWiegelmann et al 



d2005cl) and iTu et al.l d2005albl) used SUMER Doppler maps 
with potential and linear force-free models to study the outflow 
of the solar wind in coronal holes. The next natural step is to use 
more sophisticated magnetic field models, where the magnetic 
field and plasma to compute self-consistently in one model, as 
proposed here. In this paper we first present the basic equations 
used in Sect. [2] (supplemented by a couple of appendices). We 
then present the algorithm based on the equations (Sect. [3]) and 
the derivation of the non-force-free MHD equilibria we used 
for testing the code (Sect. [4]). The results are presented in Sect. 
[5] followed by a discussion and conclusions in Sect. [6] 



The MHD equilibrium equations (here without gravity for sim- 
plicity 0) are 



(VxB)xB- / u Vp = 
V B = 0, 



We have 



(1) 
(2) 
(3) 



B • Vp = 0, 

i.e. the pressure is constant along magnetic field lines. We show 
in Appendix lAl that it is possible in principle to include field- 
aligned incompressible flows into the method. We take this pos- 
sibility into account by replacing the plasma pressure p by a 
generalized pressure Yl - p + pv 2 /2 from now on. We did not, 
however, calculate any example cases with field-aligned flow, 
because one would need more information to disentangle the 
contributions of the plasma pressure and the energy density of 
the flow. 

The general form of the MHD equilibrium equations is 
given by 



(VxB)xB = VA 
V ■ B = 0, 

with the special cases 



(4) 
(5) 



A = 



po = const 
V(pop) 

v, w>n 



1 



(force-free equilibria) 
(MHS equilibria) 

(field-aligned flow, constant Ma) 



■(6) 



To solve Eqs. (01 and (0, we define the functional 



L(B, A) 



f ^1 
J B 2 



|(V x B) x B - VA| Z + w b \V ■ B| 2 d 3 x 



(7) 



where w a and Wb are positive definite weighting functions 0. 

It is obvious that Eqs. (@]i and (O are satisfied if the func- 
tional dT) reaches its minimum at L = 0. The functio nal (jT} 
generalizes the force-free approach of I Wheatland et al" ( 2000l) 
and th e magnetohydrostatic model of IWiegelmann & Inhester 
d2003l) . 

To obtain evolution equations for the magnetic field and the 
generalized plasma pressure, we take the derivative of Eq. (01 



1 While it is, in principle, possible to include gravity in the opti- 
mization principle, it is hard to find (semi)-analytic equilibria to test 
the code. For configurations that are small compared to the gravity 
scale height of some 0.1 solar radii, gravity might be neglected in 
first order, however. A consistent treatment of large-scale (some solar 
radii) structures, like helmet streamers, require not only the inclusion 
of plasma pressure and gravity, but a lso the use of spherical geome- 
try, see e.g jwiegelmann et al.1 J 1998h . but such computations are well 
beyond the scope of this paper. 

2 The functions w a an d w/, can e.g. be us ed to deal with unknown 
boundary conditions (see lWiegelmannf2004l) . To test the method, it is 
sufficient to use w a = Wb = 1 . 
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with respect to an artificial parameter t, assuming that both B 
and A depend on t : 



1 dL 

2 dt 



I 
I 



<9B 



F - V ■ fi a 



dA 



d x 



v dt " dt 

dB ~ JA 
— G + £l a n — d x, 



dt 



dt 



(8) 



where F and G are defined in AppendixlBl 

If B and A are kept fixed on the boundary of the computa- 
tional box, the surface integral vanishes and we can minimize 
L by the solving the equations 



dB 

dA 

dt 



— = -vv-a a 



(9) 



(10) 



iteratively, with positive constants p and v (see Appendix iBl 
for the mathematical derivation). The form of Eqs. © and ( fTOl ) 
ensures that L decreases monotonically during the computation 
of B and A. 

3. Algorithm 

To compute nonlinear, selfconsistent 3D-MHD equilibria, we 
use the following steps: 

1. Compute a potential field. This can be done from the B : 
component of the vector magnetogram alone. 

2. Distribute the plasma (or, say, the generalized plasma pres- 
sure A) along the potential magnetic field by solving 



B ■ VA = 



(ID 



with an upwind method (used here) or a magnetic field line 
tracer 0. 

3. Substitute the boundary values of the computational box 
from the observed vector magnetogram. The interior of the 
computational box remains filled with a potential field and 
a corresponding plasma distribution. 

4. Iterate for B and A by Eqs. (O and ( [Tol l. The continuous 
form of these equations ensures that L is monotonically de- 
creasing. This is also ensured in the discretized form if the 
iteration step dt is small enough. The code automatically 
controls the optimal iteration step. If L(t + dt) > L(t), the 
step is refused and repeated with dt reduced by a factor of 
two. We increase dt by a factor 1.01 after each successful 
iteration step to allow dt to be as large as possible with re- 
spect to the stability condition. 

5. The iteration stops when L reaches its minimum. In practise 
we stop the iteration when %/L < 1CT 6 for 100 consecutive 
steps. 

6. As result we get the magnetic field B and generalized 
plasma pressure A, which fulfill the MHD equations and 



3 By multiplying Eq. ((6) with B, we get B ■ VA = 0, which implies 
that the generalized pressure A is constant on magnetic field lines. 
This kind of equation B ■ Vc = is also used in Grad-Rubin like 
extrapolation codes to distribute the force-free parameter a along the 
field lines in space. 



is consistent with the observed boundary conditions. In the 
final step one has to disentangle the generalized plasma 
pressure A with respect to the plasma pressure p and 
the flow velocity v. For MHS equilibria, this is trivial: 
p = -f^; for equilibria with flow one needs further obser- 
vations/assumptions regarding the plasma flow, e.g. from 
SOHO/SUMER. 



4. Semi-analytic test equilibrium 

We test our method with the h elp of a sem i analy tic equi- 
librium similar to the force-free iLow & Loul ([1990) solution 



(henceforth LL). We solve the Grad-Shafranov equation for 
axis-symmetric force-free fields in spherical coordinates r, 8, 



d 2 A 1 - p 2 d 2 A 



dA , dbj, 



dr2 r „ dM ^^d-p^ +b ^. (12) 



dA v dA 



LL looked for self-similar (in r) solutions by choosing b$(a) = 
cA q . We generalize this approach by also choosing A(A) = kA s , 
where the power s has to be chosen such that self-similar solu- 
tions are possible. 

Following LL we assume 



A(r,p) = 



(13) 



Substituting this into Eq. (fT2l) and using the above expressions 
for A and b^, we get 



n(n + 1)- 



1 - p 1 d 2 P 
~diJ 



Jl + 2 yJl + 2 



p sk(l -p ) 



^(5-15-2 



+ qc- 



plq-\ 
■ndq-l) ' 



(14) 



We now determine q and s such that all powers of r are 
equal, obtaining the equations 



n + 2 = n(s-\)-2 
n + 2 — n(2q — 1). 

Solving Eq. ( fT5l ) gives 

5=2+-, 

n 

whereas Eq. ( TTol i gives the same result as LL 

q = 1 + -. 
n 

The equation for P{p) is then given by 
d 2 P I 2 



(15) 
(16) 



(17) 



(18) 



(1 - p)—j + n(n + 1)P + 2k\l + -1(1 ~p)P 



dp 



,2s. p l+4/;i 



1 



+ fl / |l + -|P 1+// "=0. 



(19) 



This equation is nonlinear and has to be solved numerically. 

For our test equilibrium, we use a similar parameter set to 
the one in LL for the force-free case (<B = | , / = 0.3). In Eq. 
dT9l > we choose (as LL) a 2 — 0.425. Both O and / have the same 
meaning as in LL here. 
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Original Reconstruction 




Fig. 1. Left panel: original MHD-equilibrium; right panel: reconstruction. The colour coding shows the line-of sight magnetic 
field on the photosphere. 



Pressure (O rig . 




Pressure (Rec.) 



1 20 30 40 50 60 
Log ( P ress u re) (Orlg. 




60 
50 
40 
30 
20 

1 

O 



1 20 30 40 50 60 
Log ( P ress u re) (Rec.) 




1 20 30 40 50 60 



1 20 30 40 50 60 



Fig. 2. We show the line-of-sight integration of A along the z-axis (Top: linear scaling, bottom logarithmic scaling). This quantity 
seems to be important because coronal images have a line-of-sight integrated character as a consequence of the optical thin 
coronal line emission. The left hand panels correspond to the original solution and the right hand panels to our reconstruction. 
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O 1x10^x1O%x10 6 4x1O%x10%x1O 6 
teration steps 

Fig. 3. Evolution of the functional L(B, N) during the iteration 
for a box of 64 3 . The dotted line corresponds to the discreti- 
sation error of the analytic solution [L = 0.16 or Logio(L) = 
-0.80]. 



The difference of our solution to LL is, that we have the 
additional term 2k(l + 2 J (1 - /i 2 )P l+4/ " which corresponds to 
the (generalised) plasma pressure, (k — corresponds to LL.) 



2.097 for 



For the present paper we choose k — 10 and j- 
At = -1. 

5. Results 

In Fig. Q] we compare some magnetic field lines of the original 
solution (left panel) with the result of our reconstruction (right 
panel). The reconstructed solution obviously agrees well with 
the original, and an inspection by eye only shows hardly visible 
differences for some loops. Figure [2] contains images produced 
from the original (left panels) and reconstructed (right panels) 
plasma pressure. The images have been produced by a line- 
of-sight integration along the z-axis to mimic the optically thin 
coronal plasma. The top panels use a linear scaling and the bot- 
tom panels a logarithmic one. While the overall structure of the 
plasma emission agrees in the original and the reconstruction, 
there are small deviations (visible in the logarithmic images) in 
darker (weak magnetic field) regions, in particular at the right 
side of the image. To evaluate the quality of reconstruction we 
used several figures of merit. The functional L as defined in Eq. 
(O provides automatically a measure of how well the MHD 
equilibria and the solenoidal condition are satisfied. 

The evolution of L during the optimization process is 
shown in Fig. [3] as a function of the iteration steps. For the 
64 3 box, L has decreased below the discretization error of the 
analytic solution during the iteration. The number of itera- 
tion steps (and thus the the computing time) for convergence 
is, however, about 1000 times greater than for a correspond- 
ing force-free optimization. The reason for this behavior is be- 
cause the semi-analytic solution has a huge difference in mag- 
nitude of the pressure throughout the computational box with 
Pmax = 1-9133 ■ 10 4 and p min = 1.0239 ■ 10 -21 . The code obvi- 



ously has problems dealing with the pressure differences during 
the iteration. Test runs with prescribed pressure profile (only 
iterating for B) converge at th e same rate as a corr esponding 
force-free calculation (see e.g, Schriiver et alj|2006l where the 
ma gnetic pressure varie s between 4.616- 10 4 and 2.87- 10~ 2 .) 



Schriiver et al. ll200d) developed several figures of merit 



to quantify the difference between two vector fields B (semi- 
analytic test field) and b (reconstructed fields). The figures have 
been used to evaluate the quality of six different non-linear 
force-free extrapolation codes, by comparing the result with the 
L L solution. We use the same definitions as given in section 4 
of ISchriiver et al. I (120061) for evaluating the quality of the non 
force-free magnetic fields here: 



- Vector correlation 



( v /2 

C vec = 2Bi-bi/ Z |Bi|2 Z |bil I ; 
Cauchy-Schwarz inequality 

JV 4-HBiHbif 

where N is the number of vectors in the field; 
Normalized vector error 

£N = 2 lbl ~ Bll/ Z |Bil; 



(20) 



(21) 



Mean vector error 

Ibi-Bil 

N - y 



—z 



|Bi| 



(22) 



(23) 



- Magnetic energy of the reconstructed field normalized with 
the energy of the input field 



(24) 



The two vector fields agree perfectly if the figures of merit 
(Cyec, Ccs, e ) aie unity and if (En, Em) are zero. 

We also compare how closely the generalized plasma pres- 
sure A of our reconstruction agrees with the original. To do 
so we compute the linear Pearson correlation coefficient, both 
for the 3D plasma pressure CorrK^ and the line-of-sight 
integration (with respect to the z axis) CorrA2D- The latter 
value is a useful quantity because observed EUV-images (e.g. 
SOHO/EIT, TRACE) have a line-of sight integrated character 
as a consequence of the optically thin coronal plasma. 

6. Discussion and conclusions 

We have genera lized the optimization method for nonlinear 
force-free fields ( [Wheatland et al J2000I) to magnetohydrostatic 
equilibria by including a pressure gradient. Using a semi- 
analytical magnetohydrostatic equilibrium s imilar to the non- 
linear force-free equilibria by lLow & Lou dl990h for testing 
the code, we showed that the optimization method also works 
in principle if plasma pressure is included. The reconstructed 
solution agrees well with the exact solution. For application 
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Table 1. Quality of the reconstruction. 





Exact 


Rec. 


L 


0.16 


0.035 


Vector Correlation 


1 


0.999 


Cauchy Schwarz 


1 


0.992 


Normalized Vector Error 





0.130 


Mean Vector Error 





0.073 


Relative Magnetic Energy e 


1 


0.996 


CorrA 3D 


1 


0.998 


CorrK^o 


1 


0.999 



to real data, the method will have to be developed further to in- 
clude measured information about the plasma properties, which 
in the present paper we have taken from the exact solution. We 
have also shown that it is possible to include an incompressible 
field-aligned plasma flow. 

Compared to the corresponding nonlinear force-free code, 
however, the MHD-optimization method presented here (for 
B and A simultaneously) is about a factor of 10 3 slower than 
force-free computation. The reason for this seems to be the ad- 
ditional equation for updating the plasma pressure. If we pre- 
scribe the correct exact pressure and solve only for the mag- 
netic field, the convergence speed of the method is similar to 
that of the corresponding force-free case. If, on the other hand, 
we fix the magnetic field using the exact solution and iterate 
only the pressure, the convergence speed is much less and com- 
parable to the convergence speed of the combined magnetic 
field/plasma pressure iteration. The reason for this behavior by 
the method is the huge difference in the values of the plasma 
pressure throughout the computational box, combined with the 
relatively low values of the plasma ft. This means that even 
small changes in the magnetic field have to compensate for by 
much larger changes in the plasma pressure. 

Despite these practical difficulties, we believe that the 
method has a lot of potential to improve magnetic-field recon- 
struction by including more information from observations. 

Appendix A: Including field-aligned flow. 

For completeness, we give some details here of how field- 
aligned incompressible flows could in principle be included 
in the optimization method. Although flows seem to be com- 
mon in the corona, the inertial force scales with the square 
of the Alfven Mach number, which usually means that only 
flows with Alfven Mach numbers close to 1 have a noticeable 



effect upon the equilibrium structure (e.g. iPetrie & Neukirch 
1999h . We also notice that to use the method including flow 
one would need additional information, to be able to disentan- 
gle the effects of the plasma pressure and plasma flow. The 
assumption of an incompressible field-aligned flow is an ide- 
alization, because any real flow that approaches the base of the 
coronal loops, which has small pressure scale heights, will be 
very compressible, even when the flow is slow. 

We restrict our treatment here to incompressible field- 
aligned flows with constant Alfven Mach number M A . For in- 
compressible fie ld-aligned flows, a general tra nsformation the- 
ory exists (e.g. iGebhardt & Kiesslind 1 19921) that allows the 



mapping of equilibria with field-aligned incompressible flow 
(with M A < 1) onto static MHD equilib ria. For further details 
we refer to lGebhardt & Kiesslind dl992l) . 

Neglecting gravity for simplicity, the stationary incom- 
pressible MHD equations are 



p(v ■ V)v 

V B 

V ■ (pv) = 
V-v = 
V x (v x B) = 



1 

/'o 




■(VxB)xB- Vp 



(A.l) 

(A.2) 
(A.3) 
(A.4) 
(A.5) 



Equations (IA.3b and iAAi imply that the plasma density is con- 
stant along magnetic field lines: 

B ■ Vp = 0. 



Equation dA.5t is identically satisfied for field-aligned flow (v || 
B, implying a vanishing electric field). The plasma velocity can 
then be written as 

v = M A v A 

where M A is the Alfven Mach number and va the Alfven ve- 
locity, defined by 



B 



va = 



#4)P 



Rewriting the force balance equation with the vector identity 

1 9 

(v ■ V)v = -Vv 2 + (V x v) x v, 
it takes the form 

p(Vxv)xv = — (VxB)xB-vL + — ), (A.6) 
W> \ 2 ) 

and for a constant Alfven Mach number M A we immediately 
get 



(VxB)xB = V 
V B = 



jU()II 

\-M\ 



(A.7) 
(A.8) 



where II = p + ^- is the generalized pressure (plasma pres- 
sure and dynamic pressure). Equation ( IA.7b has (at least for 
M A < 1) a structure that is similar to the magnetohydrostatic 
equilibrium equation. This is a natural property of MHD equi- 
libria with an incompressible stationary plasma flow, which 
can be derived from MHS -equi libria by a suitable tran s forma - 
tion as shown, for example, in IGebhardt & Kiesslind (119921) . 
Equations dA.7b and dA.8b can then, in principle, be solved 
using the method presented in Sect. 12 but more information 
would be needed to obtain the plasma pressure p and the den- 
sity p and velocity v separately. 
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Appendix B: Mathematical details 

With 



(V x B) x B - V 



n a = b 

Q b = B 2 [(V B) B] 
the functional (01 reads 



( w a B 2 Q? a + w b B 2 Q. 2 b cPx. 
Jv 



l-M 2 ) 



(B.l) 
(B.2) 

(B.3) 



We minimize equation (IB. 31 with respect to an iteration param- 
eter t and obtain an iteration equation for the magnetic field 

1 dL r <9B - , r SB ~ , 



F = F a + F b 
G = G a + Gb 

F a = W a F a + (fi a X B) X Vw a 

F~ b = w b F b + (£2 b ■ B) Vw h 
G a = w« G a 
Gb = Wb Gb 

F a = V X (Q a X B) - fi a x (V x B) + Cl 2 a B 
F b = V(Q b • B) - fi b (V • B) + Q 2 B 



G a = n x (fi a x B) 
G b = -n(n b B) 



(B.5) 
(B.6) 

(B.7) 
(B.8) 
(B.9) 
(B.10) 

(B.ll) 
(B.12) 
(B.13) 

(B.14) 
(B.15) 



and h is the inward unit vector on the surface S . The surface 
integral in ( IB .41 ) vanishes if the magnetic field is described on 
the boundaries of a computational box. 
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